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Contact Stress Analysis of Spiral 
Bevel Gears Using Finite 
Element Analysis 

A procedure is presented for performing three-dimensional stress analysis of spiral bevel 
gears in mesh using the finite element method . The procedure involves generating a 
finite element model by solving equations that identify tooth surface coordinates . 
Coordinate transformations are used to orientate the gear and pinion for gear meshing. 
Contact boundary conditions are simulated with gap elements. A solution technique for 
correct orientation of the gap elements is given . Example models and results are 
presented. 


Introduction 

Spiral bevel gears are used to transmit power between 
intersecting shafts. One such application is in helicopter 
transmission systems. In this critical application, the gears 
operate at relatively high rotational speed and transmit sub- 
stantial power (i.e., 1500 HP at 21,000 rpm). 

Prior research has focused on spiral bevel gear geometry 
[1-6] to reduce vibration and kinematic error, improve manu- 
facturability and improve inspection. Stress analysis is an- 
other important area of ongoing research. Accurate predic- 
tion of contact stresses and tooth root/fillet stresses are 
important to increase reliability and reduce weight. Tooth 
flexibility will shift the contact zone and path, alter and 
change the contact area and stresses. Finite element analysis 
of gears in mesh, including tooth flexibility, will model con- 
tact stresses better than a Hertzian calculation. 

Much effort has focused on predicting stresses in gears 
with the finite element method. Most of this work has in- 
volved parallel axis gears with two dimensional models. Only 
a few researchers have investigated finite element analysis of 
spiral bevel gears [7, 8]. 

Finite element analysis has been done on a single spiral 
bevel tooth using an assumed contact stress distribution [9]. 
The research reported here will utilize the numerical solu- 
tion for spiral bevel surface geometry to study gear meshing. 
Pinion tooth and gear tooth surfaces will be developed based 
on the gear manufacturing kinematics. The individual teeth 
are then rotated in space to create a multi-tooth model (4 
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gear and 3 pinion teeth). The tooth pair contact zones are 
modeled with gap elements. The model development proce- 
dure and finite element results are presented. 


Equations for Tooth Surface Coordinates 

The system of equations, briefly summarized here, re- 
quired to define the coordinates of a face-milled spiral bevel 
gear surface were developed by Handschuh and Litvin [9]. 

A conical cutting head, attached to a rotating cradle, 
swings through the work piece. Parameters U and 0 locate a 
point on the cutting head in coordinate system S c attached to 
the cutting head as shown in Fig. 1 and described by the 
following equations. 

r cot ip - U cos tp 
U sin ip sin 9 
U sin ip cos 0 
1 

The roll angle of the cradle <p c is used to locate the rotating 
cradle with respect to the Fixed machine coordinate system 
S m . Parameters U, 6 and <p c , along with various machine tool 
settings can be used to completely define the location of a 
point on the cutting head in space. 

Since the kinematic motion of cutting a gear is equivalent 
to the cutting head meshing with a simulated crown gear, an 
equation of meshing can be written in terms of a point on the 
cutting head (i.e., in terms of U, 6 and <b c ). The equation of 
meshing for straight-sided cutters with a constant ratio of roll 
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relative to the cutting head and cradle coordinate systems (S c 
and S s ) respectively. These parameters can be transformed 
through a series of coordinate transformations to a coordi- 
nate system attached to the work piece. Or U, 6 , <p c can be 
mapped into X w , Y w > Z w in coordinate system S K , attached to 
the work piece. These transformations, used in conjunction 
with two other geometric requirements, give the two addi- 
tional equations. 

The correct V, 0 and <p c that solves the equation of 
meshing, must also, upon transformation to the work piece 
coordinate system 5 H „ result in an axial coordinate Z H , that 
matches with a preselected axial position Z. (See Fig. 2) 

Z K . - Z = 0 (4) 


between the cutter and work piece is given in [1, 5, 6] as: 

(U - r cot tp cos ip ) cos y sin r + 5[(/n w - sin y) 

Xcos ip sin 8 + cos y sin y sin (<7 - <p c )] 

±£ m ( cos y sin ip + sin y cos ip cos r) 

-L m sin y cos ip sin r = 0 (2) 

The upper and lower signs are for left and right hand gears, 
respectively. The following machine tool settings are defined 
[5, 6, 9]. 

i p - cutting blade angle 
r = (0 + q ± <p c ) 
q - cradle angle 
y = root angle of work piece 
E m = machining offset 

L m = vector sum of change of machine center to back and 
the sliding base 

m cw = <p c /(p w , the relationship between the cradle and work 
piece for a constant ratio of roll 
U = generating cone surface coordinate 
S = radial location of cutting head in coordinate system 

r = radius of generating cone surface 
This is equivalent to: 

f } (U, 6 , 4> c ) = 0 (3) 

Because there are 3 unknowns, U , 6 , and <p c , three equations 
must be developed to solve for the surface coordinates of a 
spiral bevel gear. The 3 parameters U, 0 and <p c are defined 


This equation along with the correct coordinate transforma- 
tions [see Eq. (11)] result in a second equation of the form: 

/ 2 (t/, 0, <p c ) = 0 (5) 

A similar requirements for the radial location of a point on 
the work piece results in the following: 

^ - {*:■ + Y;)' r ~ = o (6) 

The appropriate coordinate transformations (see Eq. (11)) 
convert Eq. (6) into a function of (/, 0, and <p c . 

/ 3 (U, 0, <p c ) = 0 (7) 

Equations (3), (5) and (7) form the system of nonlinear 
equations necessary to define a point on the tooth surface. 


Solution Technique 

An initial guess f/°, 0°, <p° is used to start iterative 
solution procedures. Newton’s method is used to determine 
subsequent values of the updated vector (£/\ 0 A , <£*)' [10]. 

Where the vector Y is the solution of: 

The 3x3 matrix in the preceding equation is the Jacobian 
matrix and must be inverted each iteration to solve for the Y 
vector. The equation of meshing, function is numerically 
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differentiated directly to find the terms for the 


( 8 ) 


y*- 1 ' 

y 2 *-‘ 

\rk — 1 


= (9) 

Jacobian matrix. Function f 2 and / 3 cannot be directly 
differentiated with respect to U, 6 and <f> c . After each itera- 
tion U k ~\ 6 k ~\ 4>c 1 (in the cutting head coordinate 
system Sc) are transformed into the work piece coordinate 
system, S K „ with the series of coordinate transformations as 
given in Eq. (11). 

{ r cot ip — U cos ip\ 

V sin (A sin <9 l (10) 
U sin if/ cos 6 } 




y*-l 

Y k ~ x 


df t (U k ~') 

w*- 1 ) 


'A 

dU 

dd 

d(f> c 



df 2 (6 k -') 


A 

dU 

duQ 

d<t> € 


df,{u k -') 


* flirt- 

A 

dU 

de 

84> c 



where: 


[M„ c ] = [M K J(4> c )][M ap ][M pm ][M ms f(<t> c )][M sc ) 

( 11 ) 

Each matrix [M] above represents a transformation from one 
coordinate system to another. (See Appendix I for the spe- 
cific matrices.) 

Functions f 2 and / 3 are evaluated by starting with an 
initial U k , 9 k and <}> k , performing the transformation in Eq. 
(11) and evaluating Eqs. (4) and (6). The numerical differen- 
tiation of f 2 and / 3 is performed by transforming U k + inc, 
6 k + inc, <f> k + inc (where inc is a small increment appropri- 


ate for numerical differentiation into X w + inc, Y w + inc, 
Z w + inc. Equations (4) and (6) are then used to evaluate the 
numerical differentiation. Function /,, / 2 , / 3 and the partial 
derivatives of /„ / 2 , / 3 required to the Jacobian matrix are 
updated each iteration. The iteration continues until the Y 
vector is less than a predetermined tolerance. This completes 
the solution techniques for a single point on the spiral bevel 
gear surface. 

The four comers of the active profile are identified from 
the tooth geometry plane as shown in Fig. 2. Point 1 on the 
surface is chosen to be the lowest point of the active profile 
on the toe end. The initial guess to start the procedure has to 
be sufficiently close to the correct solution for convergence to 
occur. The solution for the first point proves to be an 
adequate initial guess for any subsequent points on the 
surface. 

Subsequent interior surface points are found by increment- 
ing r = (X 2 + Y 2 )^ and Z. By adjusting the increments 
used, a surface mesh of any density can be calculated. The 
process is repeated four times for each of four surfaces; gear 
convex, gear concave, pinion convex and pinion concave. 
Software was written to solve for the tooth surface coordi- 
nates. Additional software converted the surface coordinates 
into input commands for the 3D solid modeler PATRAN 
[11]. PATRAN was used to create the nodes and elements 
for the FE code NASTRAN [12]. 

Since all four surfaces are generated independently, addi- 
tional matrix transformations are required to obtain the 
appropriate orientation for meshing. The proper convex and 
concave surface orientation (for both the gear and pinion) is 
found by fixing the concave surface and rotating the convex 
surface until the correct tooth thickness is obtained. The 
correct angle of rotation is obtained by matching the tooth 
top land thickness with the desired value. 


Gear and Pinion Orientations Required For Meshing 

After generating the pinion and gear surface as described 
above, the pinion cone and gear cone apex will meet at the 
same point as shown in Fig. 3. This point is the origin of the 
fixed coordinate system attached to the work piece being 
generated. To place the gear and pinion in mesh with each 
other rotations described in the following example are re- 
quired: 

(1) The 19th tooth is selected for meshing. This corre- 
sponds to rotating the generated tooth 190 deg. CW about 
the Z w axis. (For this example, N , *= 36 total teeth on the 



Fig. 3 Orientation of gear and pinion based on solution of the system of 
equations, and after rotations required for mesh 
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Fig. 4 Complete 3D model of gear and pinion in mesh 


gear.) For the general case the gear tooth is rotated 360/7V, + 
180 deg CW about its axis of rotation (Z H ,). This corresponds 
to selecting the ith tooth of the gear to be in mesh where 

i = s t /2 + 1. 

(2) The pinion is rotated by 90 deg. CCW about the Y 
axis. Note: rotation (1) corresponds to selecting a different 
tooth on the gear to be in mesh; however, rotation (2) 
corresponds to physically rotating the entire pinion until it 
meshes with the gear. 

(3) Because the four surfaces are derived independently, 
their orientation is random with respect to meshing. The 
physical location of the gear and pinion after rotations (1) 
and (2) correspond to the gear and pinion in mesh with 
severe interference. To correct the interference the pinion is 
rotated CW about its axis of rotation (Z M ) until surface 
contact occurs. For this example the rotation was 3.56 deg. 
Figure 4 shows an example of a simulated gear pair meshing. 
The generated pinion tooth was copied and rotated 12 times 
and the generated gear tooth was copied and rotated 36 
times. 

Contact Simulation 

The tooth pair mesh contact point can be located by a 
method described by Litvin [6] or by a search technique. 
Pairs of finite element node points (one on each tooth 
surface) are evaluated until the pair with smallest separation 
distance is obtained. (A finer finite element mesh would yield 
greater resolution.) Once the contact point is established, a 
vector normal to the surface at the contact point is calcu- 
lated. 

The intersection of the normal vector on the pinion at the 
contact point with the gear surface identifies the second 
point required to define the gap element. Additional gap 
elements are obtained by taking additional vectors from 
other pinion surface finite element nodal points (parallel to 
the contact point normal vector), and calculating where they 
intersect the mating gear surface. Finite element nodal points 
on the gear surface are located to the intersection points of 
the normal vectors and the gear tooth surface. 

The vector normal to the cutting surface is given in [5] as: 

( sin ) 

{O = | cos iJj c sin t l (12) 

I COS COS T I 

This vector is written relative to S r , a coordinate system fixed 
to the cutting head. To obtain the vector n M relative to the 
coordinate system fixed to the work piece, the following 
transformation must be performed. 

{-,} = [L~c){ nJ (13) 

Where [L wc ] is found by removing the 4th row and column 
from [M wc ]. 



Fig. 5 Distorted gear after connecting gap elements 


Each gap element is connected between a node on the 
pinion surface, hereafter designated node 1, and the corre- 
sponding intersection node point on the gear surface, desig- 
nated node 2. The intersection point on the gear surface is 
found as follows: 

Consider node 1, a point on the pinion with coordinates 
X p> Y , Z p in S w . Let Q x , Q yy Q z be any point in space such 
that ( Q x , Q y , Q z )‘ = ( X p , Y py Z p ) 1 + b(n x , n y , n z )‘ where b 
is a scale factor, and n x , n y and n z are the components of 
the normal vector at the contact point. The intersection of 
the normal vector from node 1 with the mating gear surface 
defines node 2 and has to satisfy the following three equa- 
tions: 

Qx = G x “ X P + bn x 

Q y = Gy = Y p + bn y 

Q: = = Z p + bn z (14) 

Where G x , G v , G z is a point on the gear surface. A point on 
the gear surface must also satisfy 


l G A 

1 

(r cos - U cos <// f \ 

Cvl 

GJ 

= [W w ] 

U sin i// c sin 0 1 

U sin i Ij c cos 0 I 

1 J 


1 1 ] 


Equations (14) and (15) lead to a system of 3 nonlinear 
algebraic equations. These three equations, along with the 
equation of meshing for the gear surface, provide a system of 
4 equations and 4 unknowns (w, 0, <t> c , b). These equations 
are solved with Newton’s method described earlier. The 
intersection of the normal from node 1 on the pinion with 
the gear surface is now obtained. 

This procedure is used to locate the intersection of nor- 
mals from all points on the pinion surface (in the contact 
zone) with the gear surface. The gear tooth surface is 
remeshed utilizing the intersection points as shown in Fig. 5. 
Gap elements are connected between corresponding nodal 
points on the pinion and the intersection points on the gear 
surface. 

Finite Element Model 

An example analysis was performed using gears from the 
NASA Lewis Spiral Bevel Gear Test Facility. In this case, the 
left hand pinion mates with the right hand gear. Counter 
clockwise rotation of the pinion results in the concave surface 
on the pinion mating with the convex gear surface. 

The design data for the pinion and gear are given in Table 
1. The design data are used with methods given in [6] to 
determine the machine tool settings for the straight sided 
cutter data given in Table 2. 

The finite element gear pair model, shown in Fig. 6, 
contains 4 gear teeth and three pinion teeth. The model had 
10,101 nodes (30,303 degrees of freedom) and 7596 eight 
noded 3D brick elements. The analysis was done on a gen- 
eral purpose FE code [12]. 
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The three pinion teeth are fixed in space with zero dis- 
placements and rotations. This is done by setting jc, y and z 
displacements equal to zero on the four comer nodes of each 
rim section. The gear is constrained to rotation about its axis 
of rotation. The gear is loaded with a torque of 9450 in lbs on 
the gear by applying 4725 lb force located 2.0976 from the 
gear axis of rotation. 

At the orientation chosen between the pinion and gear 
two pairs of teeth were in contact. One pair had contact near 
the middle region of the tooth and another pair had contact 
near the toe (i.e., about to go out of mesh). Initially the 
model started with a total of twenty one gap elements. For 
the tooth that is approximately midway through mesh. Fifteen 
gap elements were used. For the tooth about to leave mesh, 
four gap elements were used. The analysis starts with one 
gap elements closed in each contact zone. Within the finite 
element code an iterative process is used to determine how 
many gap elements must close to reach static equilibrium. 
The solution iterated four times before reaching equilibrium 
with four gap element closed in the main contact region and 
one closed in the edge contact region. Stress contours for the 
pinion tooth with contact are shown in Fig. 7. 

The average nodal minimum principal stresses in the main 
contact zone average -204,000 psi with a maximum of 
-299,900 psi and a minimum of -103,574 psi. The corre- 
sponding elemental stresses average - 103,500 psi with a 
maximum and minimum of -123,900 psi and -79,500 psi, 
respectively. The nodal stresses are higher because of load 
concentration from the gap elements. These stress ranges 
bracket the estimated hertzian stresses for the gear set under 
the same load conditions. Contact with only 4 gap elements, 
along with large stress gradient among adjacent nodes indi- 
cate the need for a finer finite element mesh for improved 
stress prediction. 


Table 1 Pinion and gear design data 



Fig. 6 Seven tooth model used for finite element analysis of mesh 



Fig. 7 Stress contours in pinion tooth 



PINION 

GEAR 

Number of teeth pinion 

12 

36 

Dedendum angle, deg 

1.5666 

3.8833 

Addendum angle, deg 

3.8833 

1.5666 

Pitch angle, deg 

18.4333 

71.5666 

Shaft angle, deg 

90.0 

90.0 

Mean spiral angle, deg 

35.0 

35.0 

Face width, mm fin) 

25.4(1.0) 

25.4 (1.0) 

Mean cone 
distance, mm (in) 

81.05(3.191) 

81.05 (3.191) 

Inside radius of gear 
blank, mm fm) 

5.3 (0.6094) 

3.0 (.3449) 

Top land thickness, mm fm) 

2.032 (0.080) 

2.489 (.098) 

Clearance, mm On) 

0.762 (0.030) 

0.92964 (0.0366) 


Conclusions 

A multi tooth finite element model was used to perform 
three-dimensional contact analysis of spiral bevel gears in 
mesh. Four gear teeth and three pinion teeth are generated 
by solving the equations, based on gear manufacturing kine- 
matics, that identify tooth surface coordinates. The gear and 
pinion are orientated for meshing with coordinate transfor- 
mations. Software was written to solve for the tooth surface 
coordinates, and create input commands for a 3D solid 
modeler. The solid modeler created the nodes and elements 
used as input for a general purpose finite element code. 

Surface stresses are evaluated with gap elements. The gear 
surface is remeshed with nodal points relocated to identify 


Table 2 Generation machine settings 


PINION 


GEAR 


Radius of cutter, r, in 
Blade angle, deg 
Vector sum, 

Machine offset, E., in 
Cradle to cutter distance, s, in 
Cradle angle q, deg 
Ratio of roll, 

Initial cutter length, u, in 
Initial cutter orientation, 6 , deg 
Initial cradle orientation, deg 


Concave 

2.9656 

161.9543 

0.03849 

0.15457 

2.94780 

63.94 

0.30838 

9.59703 

126.83544 

-0.85813 


Convex 

3.0713 
24.33741 
-0.05181 
-0.17426 
2.80104 
53.926 
0.32204 
7.42534 
124 43689 
-11.38663 


Concave 

3.0325 

58.0 

0.0 

0.0 

2.85995 

59.23420 

0.95086 

8.12602 

233.9899 

-0.35063 


Convex 

2.9675 

22.0 

0.0 

0.0 

2.85995 

59.23420 

0.95086 

7.89156 

234.9545 

-12.3384 


Journal of Mechanical Design 


JUNE 1995, Vol. 117/239 




the correct gap element orientation. Initial FEA stress results 
compare favorably with calculated hertzian contact stresses. 
However; large stress gradients between adjacent nodes in 
the contact zones indicate a need for greater finite element 
mesh refinement. 
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APPENDIX I 

Coordinate transformations involving both rotation and 
translation require mixed matrix operations of multiplication 
and addition. Matrix representation of coordinate transfor- 
mations will need only multiplication of matrices if position 
vectors are represented by homogeneous coordinates [1]. The 
following 4X4 matrices are required to transform the ho- 


mogeneous coordinates of a point on the cutting head to a 
point on the work piece. 

Matrix [ M sc ] transforms coordinate system S c , attached to 
the cutting head, into system S s , rigidly connected to the 
cradle. 


[M sc ] = 


1 

0 

0 

0 


0 

cos q 
±sin q 
0 


0 

Tsin q 
cos q 
0 


0 

Tsin q 
s cos q 
1 


(16) 


Matrix [M ms ] transforms coordinate system S s into system 
S m attached to the frame. 


[M ms ] 


10 0 0 
0 cos 4> c Tsin <t> c 0 

0 ±sin <f> c cos 4> c 0 

0 0 0 1 


(17) 


Matrix [M ap ] transforms coordinate system S m into system 
S p which orientates the apex of the gear being manufactured. 




cos 8 

0 

-sin 8 

-L m sin 8 

0 

1 

0 

±E m 

sin 8 

0 

cos 8 

L m cos 6 

0 

0 

0 

1 


Matrix [M ap ] transforms coordinate system S p into system S a 
and locates the apex of the gear being manufactured with 
respect to S m 




COS fl 

0 

sin /x 

0 

+ sin 4> w 

cos <t> w 

0 

0 

-sin /x 

0 

cos fx 

0 

0 

0 

0 

1 


Matrix [M wa \ transforms coordinate Sa into system Sw at- 
tached to the work piece. 


[ALJ 


cos «f> w ±sin 0 0 

•f sin <f> w cos tj> w 0 0 

0 0 10 

0 0 0 1 


( 20 ) 
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